Numerical Simulation of Two Dimensional Electron Transport in a Circularly Symmetric 
Cylindrical Nanostructure using Wigner Function Methods 
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We have constructed a lattice Wigner- Weyl code that generalizes the Buot- Jensen algorithm to the calculation 
of electron transport in two-dimensional circular-cylindrically symmetric structures, where the Wigner function 
equation is solved self-consistently with the Poisson equation. Almost all of the numerical simulations to date 
have dealt with the restriction of the problem to one dimensional transport. In real devices, electrons are not 
confined to a single dimension and the coulombic potential is fully present and felt in three dimensions. We 
show the derivation of the 2D equation in cylindrical coordinates as well as approximations employed in the 
calculation of the four-dimensional convolution integral of the Wigner function and the potential. We work under 
the assumption that longitudinal transport is more dominant than radial transport and employ parallel processing 
techniques. The total transport is calculated in two steps: (1) transport the particles in the longitudinal direction 
in each shell separately, then (2) each shell exchanges particles with its nearest neighbor. Most of this work is 
concerned with the former step: A ID space and 2D momentum transport problem. Time evolution simulations 
based on these method are presented for three different cases. Each case lead to numerical results consistent 
with expectations. Discussions of future improvements are discussed. 



I. INTRODUCTION 

Electron transport in a resonant tunneling structure (^TS) 
has been studied in detail over the past decade lll|3,|3l@,Qlal- 

Almost all of the numerical simulations have dealt with the 
restriction of the problem to one dimensional transport. Even 
though much progress has been made using the ID theory, in 
real devices electrons are not confined to transport in a single 
dimension and the coulombic potential is fully present and felt 
in three dimensions. Here we present a method for numer- 
ical simulation of electronic transport through a cylindrical 
device that possesses azimuthal symmetry. Figure^a) shows 
a schematic representation of such a device. 

We work under the assumption that longitudinal transport 
is more dominant than radial transport. The total transport is 
calculated in two steps: (1) transport the particles in the lon- 
gitudinal direction in each shell separately, then (2) each shell 
exchanges particles with its nearest neighbor. During a given 
time step the particles are advanced longitudinally through the 
device, as in a ID problem, but with the inclusion of radial 
momentum. This changes the form of the potential and in- 
teraction terms of the familiar ID Wigner function (transport) 
equation (WFE). Since the latter step is computationally sim- 
ple, most of this work is concerned with the former step: AID 
space and 2D momentum transport problem (lxH-2k). In this 
paper, the WFE is solved self-consistently with the Poisson 
equation. 

In order to perform the numerical simulation, parallel pro- 
gramming techniques are used. A simplest way to attack this 
problem to slice up the device into cylindrically concentric 
shells as shown in figure[nb). The two above steps now be- 
come: (1) each processor (shell) calculates the (lxH-2k) trans- 



(a) Side View (b) Top View 

Figure 1: Cylindrical RTD 



port problem, then (2) each processor exchanges particle in- 
formation with its nearest neighbor. 

This papers is organized as follows: The first three sections 
detail the derivations: Section |II| the 2D WFE in cylindrical 
coordinates (assuming azimuthal symmetry), and sections Hill 
and II VI the discretization. Section|V]discusses different meth- 
ods of solving the lxH-2k problem regarding the limitations of 
today's computational resources. This includes are-derivation 
of the potential term in the WFE and a splitting of the po- 
tential into static (conduction band edge) and changing (self- 
consistent) potentials. Finally, our concluding remarks are in 
section rV III 



II. DERIVATION 

The 3D form of the Wigner function equation (WFE), can 
be written as (scattering will be added in later) 
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The derivation of the WFE was in the same vein as Frensleyfl 
(but kept in full 3D form here), and, as he has noted, has been 
derived while assuming no boundaries exist. Because of this, 
the integral limits have been omitted in the above and will be 
discussed later. 

First, let us rewrite this equation in cylindrical coordinates 
by making the substitutions: q f^z,^ ; k ki-,k^,X<i) \ y 
p,^,9. It is important to note that while q represents the spa- 
tial (i.e.: center of mass) coordinate with respect to the origin 
on the cylindrical axis, y represents the spatial coordinate with 
respect to the origin at an arbitrar3/point within the cylinder, 
not necessarily on the axis. Figure[nb) illustrates the geome- 
try of these two variables. Equation^now becomes: 



df{r,z,<j>,ky,kr,X<t,) 



h 
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dt 

dr r d(j) ' dz 



f{r,Z,(t>,k,,kr,X,l,)^ 
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/rZJl r r rlTl 

dk'^ 1^ 141^4 I d^ I dp \p\ddy 

g-2i[{k,-K)l;+{k,cosx^,-k',cosx'^)p] 



{v{r+p,z+i;,<p + e)-v{r-p,z-^.(t>-e)}f{r,z,<t>,k',X,x',p) 

(2) 

Considering a 2D problem in cylindrical coordinates with 
azimuthal symmetry, the first thing to notice is that even 
though the potential is symmetric in {V{r,z,(l>) — > 
V{r,z)), the potential difference is dependent on angle, 
y(r±p,z±C,0±0)^y(r±p,z±C,±0). Next, integrate 
out the remaining azimuthal spatial and momentum compo- 
nents. Using 

rln rlTC 

j I fir,z,^,k„kr,X<l,)\r\d(l>\kr\dx^ = {27ty\r\\kr\f{r,z,k„kr), 

(3) 

one obtains 



df{r,zX.,kr) 



dt 
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f{r,z,k,,kr)~ 



^jdk[j dk[ fj \k',\dx'^ JdcJ dpj^^'ipide 



2n 



dXipe 



{y(z+C,/- + p,0)-y(z-C,r-p,-0)}/(r,z,<,4,4). 



(4) 



The drift term is relatively simple, but the potential term is a 
bit completed. Rewriting the potential term as 
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\p\dpJip,kr)y{z,C,r,p)^{z,r,p,ki,k'^) 



using the following definitions 



J(p,kr) = I rfX^e""*'™^^*" 



(5) 



(6) 



y{z,^,r,p) = / de{v{z+i;,r+p,e)-v{z-^,r-p,-e)} 
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(7) 



.^{z,r,p,kiX) = I d4e+'"'™^^*''/(r,z,fcU;,4), (8) 

there are some reductions in complexity. Equation|6lis just the 
definition of the Bessel function 27iJo{2krP). Equation is 
evaluated at agiven p , 6 by noting that p{9) = p cos 9, leav- 
ing equation|8|to be dealt with. By expanding the exponential 
in terms of Bessel functions of x'tj,, the integral becomes 

.^{Z,r,p,k'„k'r) = r^x'^^K X V(24p)/m('-,Z,fc;„fc;)e'("+"'''^^ 

= 27t'£j,„{2k',p)f„{r,z,k'„k',). (9) 
m 

Now there is an infinite series in m, but only the term m — 
must be counted. This is because it represents the azimuthally 
independent functions, which, as dictated by the current prob- 
lem, is the form that the Wigner distribution function, /, 
should take. 

So, finally, the complete 2D form of the WFE in az- 
imuthally independent cylindrical coordinates is 



dfir,z,k,,kr 
dt 



n 

m* 



dr dz 



f{r,z,k,,kr) 



dk[ / \p\dp^{r,p,z,kr,kz~k[).^{r,p,z,k'z). (10) 



where (notice the terms and Y have been redefined 

from what was written above) 



^{r,p,z,k,) 



\k',\dk'M2k'rP)f{r,z,kzX) 



(11) 



'^{r,p,z,k„kr) = J d^sm(2k,i;)Jo{2krPy/{zX,r,p) (12) 



y{z,^,r.p)= / de{V(z+^,r + pcos0)-V(z-^,r-pcose)} 
./() 

(13) 

/■+*:?'" r+k^;'"" rL/2 

The integral limits are / dki, / dk',., / c/^. 



_J^max 







R/2 



and / pdp for a cylindrical system of length L and radius 
J 

R, recognizing that r > always. It is worthwhile to note here 
that the since the momentum variable comes from the Fourier 
transform / t/re^'*' "", the value of k'""'^ is determined by the 
spatial length of the box. This will become important when 
the problem is discretized. 

As done by many others| 1, 4], scattering is included simply 
by the addition of a relaxation time approximation, given as 

d.f{r,z,k,,kr) 



dt 



coll 



Ufi ,, I\K\dkyfdk'j{r,zXX) . 



S\k'^\dky[dk'Joir,zXX-) 



(14) 
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where the relaxation time, T is computed from the material pa- 
rameters describing scattering due to: ionized impurities and 
longitudinal, piezoelectric, acoustic & optical phonons. 



III. DISCRETIZATION 

The WFE in discretized (matrix) form is written as ^ = 
(T + U + S)f — B, where T represents the drift (kinetic) op- 
erator, U the potential operator, S the scattering (interaction) 
operator, f the Wigner function and B the boundary conditions 
arising from the drift term. In this section, we will present the 
details of our discretization of the equations derived above. 



A. Variables, Operators and Functions 

The space and momentum variables are discretized as 



^(0 = 


i(2/-l)Az, I 




Az = L/W; 


(15) 


r{n) = 


^(2«-l)Ar, 


n= l..Nr 


, Ar = R/Nr 


(16) 


(j) (in) = 


{in — l)A(j) , m 


= l-N^ , 


A<l> = 27t/N^ 


(17) 


kzU) = 




A^z , j = 


l..Ni,_,AL^ 


Tt/AzNk, (18) 


kril) = 




)Akr,k = 


l..Nk^,AK = 


Tt/ArNkM"^) 



and the functions /, and are discretized as 



f{r,z,L,kr) - 


f{n,ij,l) 


(20) 


V{r-r',z~z') - 


V{n~n' ,i — i) 


(21) 




'2^{n,n',i,j,l) 


(22) 


f[zd.rA - 


-> i ,n,n) 


(23) 




^{n,n',i,i). 


(24) 



We note here that since p and C, are on the same grid as r and 
z, they will be denoted as r' and z', respectively. 



B. 2D Poisson Equation 

We use a Fourier transform method to solve the 2D Poisson 
equation. The Fourier transform of the charge density with 
respect to z, p7, is calculated using a fast Fourier sin trans- 
form (sinFFT). Each shell exchanges p7 to every other shell 
so that each shell knows 'p^{r), the electron density of the en- 
tire cylinder. Using a standard tridiagonal solver, the sinFFT 
of the 2D potential, (j) (r) is calculated, then transformed back 
(via a sinFFT) to the full 2D potential, ^{r,z)- 



C. Drift & Boundary Conditions Terms 

We will separate the drift (or kinetic) term into longitudinal 
and radial components and treat each one slightly differently 



[T-{]{r,zAX) = 

m* 



d d 

-1 ky 

or oz 



which is broken up into 

[Tz -f] {r,z,k^,kr) 



Tik d 

l—f{r,z,k,,k;) (26) 

m* az 



hkr d 



f{r,zAX) (27) 



First the longitudinal drift term, equation|26l is computed us- 
ing a second order "upwind/downwind" differencing scheme: 



^ =^ T 2^ [3/(x) - 4/(x ± Ax) + /(x ± 2 A x)] . 

(28) 

For k^ < 0, the upwind scheme is used, and for k^ > 0, the 
downwind scheme is used, giving 



kz^O: -^f{r,z,k,,kr)^ 
oz 

T ^ [3/(r, z, ,K)- 4/(r, z ± Az, k, , k,-) + f{r, z ± 2 A z, , fc,)] . 
2 Az 

(29) 



Which, when discretized, gives 



hAk, 



;|:^:[T..fl(,v-,;,/)-±^x 
(2; - % - 1) [3/(«, /) - Af(n, i±l, j, I) + f{n, i ± 2, j, l)\ . 



(30) 



When the second order differencing scheme gets to the bound- 
ary, it is advantageous to only have the function extend one 
unit distance into the boundary. For the upwind scheme, this 
occurs at ; = A'z — 1 and / = Nz, and / = 1 and / = 2 at 
for the downwind scheme, (the «,/ indexes will be dropped 
since there is no dependence on them). When / — l,Nz a first 
order upwind/downwind differencing scheme was employed 

(^^^^^ — T'^^^^^^^p^^ in order to preserve the continuity of 
the derivative. By saying that the distribution function past 
the boundaries have a constant value, /(/ —N^ + l,j) — f{i — 
Ojj) = ffermiU) (the two dimensional Fermi distribrution), 
these positions become constant longitudinal boundary con- 
ditions, Bz(/,y), defined as: 



fir,zA,kr) (25) 



Bz(2,;) 



Bz(A^.-l,;) 

RziNzJ) 



jffenni (j) 
^Cjf Fermi {j) 



J - 2 



Cjf Fermi {j) 
'^Cjffermiij)- 



(31) 
(32) 



(33) 
(34) 
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Where Cj = {ij-Nk.-l)- This makes the final form 

of the longitudinal drift term, [T^ • f] (/, j), as: 

Nk, 



[Tz-f](/=l,y) 
[Tz-f](2,j) 



J> 



-2Cj[2f{i=\,j)\ 

-Cj[3f{i = 2,j)-Af{i=l,j)] 

--Cj[3f{i,j)-4f{i,j)+f{i,j)] 



(35) 
(36) 
(37) 



D. Potential Term 



The potential term in eauationFTOlis written in operator form 



+ - 



1 



Tlh J-t 



-K 



[\}-t](r,z,h,K) = 

rR/l 

dki / \r'\dr'-W (r,/ ,z,k, - ki,kr)^(r.r' ,z,k'.) 



(48) 



[T,.f] (,■,;■) = Cj[3f{iJ)-4f{iJ)+f{i,j)] (38) where, 
[Tz • f] (^z - 1 , j) = Cj [3f{i = N,-l,j)- Af{i = N, , j)] (39) 
[T^-f]{N,.j) = 2Cj[2f{i = N„j)]. (40) 



Equations |^ to |33 completely define the discretized longi- 
tudinal boundary conditions, while equations 1351 to B71 com- 
pletely define the longitudinal drift for positive momenta and 
equations |38]to^]completely define the longitudinal drift for 
negative momenta. 

Next, the radial drift term, Tr ■ f, is computed using a first 
order differencing scheme since having each radial shell talk- 
ing only to its nearest neighbor is needed when this algorithm 
is parallelized. For a shell that has neighbors on both sides, 
a central differencing scheme (CDS) is used. For the inner- 
most and outermost shell, a forward/backwards differencing 
scheme (FBDS) is employed (the indexes /, j are omitted since 

there is no dependence on them) (C = ): 



(41) 



[Tr-f](n=l,/) ^ -C[/(« =!,/)-/(« = 2,0] 
[Tr-f](«,/) ^ /)-/(«-!,/)] 
[Tr-f] = ^ C[/(« = A'.,/)-/(n = A',-l,/)l. (43) 



(42) 



This form dictates that some tricks must be used at the in- 
nermost and outermost rings. For the innermost shell, it 
can be imagined that any particle possessing negative mo- 
menta (traveling inwards) will pass through the middle and 
then posses positive momenta (traveling outwards). This de- 
mands a "particle mirror" at the origin by, basically, saying 

% 

2 



y{z,Z,r,r') = j^^ d(l>' {V{z + Z,r + r'cos(j>')- (49) 
-V {z-z',r-r'cos(j)')} 

^{r,r',zA) = r''' \k[\dk'M'^k'/)f(r,z,k,X) (50) 

rL/l 

'^{r,r',z,k„kr) = / dz'sm(2k,z')Jo{2krr')y{z,z',r,r'X51) 
J 

For a given longitudinal position and radius, (z, r), one sweeps 
through disks of constant z'. For each disk of constant 
z', the potential difference contribution, V (r + r' cos — 
V {r — r' cos (j)'), between all points on the disk on (z, r) is cal- 
culated (see figure nib)). The result is that the effect of the 
potential of the entire cylinder on a given point (z, r) operates 
on the distribution function for that point. The terms 
and y each present some difficulties, but they only involve 
one integral each. It is best to go through each one separately. 

The y term, when discretized, becomes 

r{i,i',n,n') = A^ ^ {y (/+/',«+«") -y(/-/', «-«")}, 

m' = l 

(52) 

where 



that all values of < ^ j at time f, will be added to n" = n'lNT [cos(^(m')] =n'lNT [cos ( [m' - l] A(/))] , (53) 



f \\,l> -Yj at time t 
values of /(a^^, / > % 



h Af . For the outermost shell, all 
at time f, will be subtracted from 



/(.. 



/ > 



at time f + Af and, depending on the chosen 



INT[x\ being a function returning the nearest integer to the 
real value x. 

We rewrite the potential term, by introducing new func- 
tions, as 



external conditions, given values of / [Nr^l < -f-j will be 
added at each time step. This is expressed as radial boundary 

%-l)): 



conditions as (Q = (21 



[V-f\(r,z.k,,k,.) = —f I ' dk'J ' dl<!,U{r,z,k--K^K,K)f{r,Z:^X) 

(54) 



Br()l= 1,1 < 



Nk 



Nk 



Ni, 



Br{n = Nr,l< 



Nk, 



-c, 
-c, 

HQ 

HQ[+/c,-,,„(/)] 



Nk ' 

+/(n=l,/<^) 

Nk ^ 

-f(n=Nr,l>^] 



(44) 
(45) 
(46) 
(47) 



rR/l 

U(r,z,k,,kr,k',) = \k'J / dr'\r'\jo(.2ky)Jo{2krr')^(r,r',z,k.) 
Jo 



(55) 



rL/2 

3^{r,r',z,k,)= dz'sm{2Lz')y{z,z,r,r') (56) 
./ 

where Y{z,z' ,r,r') is defined as above. In discretized form. 
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these equations become 

[U-f](n, /,;,/) 
U(n.iJ,IJ') 



Y^Y^U{nJJ-fjj')f{nJj'j') (57) 



(2/'-% - l)\Y^J{n',l,l')^P{n,n',i,j) (58) 

n' 

(2«' - 1 ) 1 7o ^ (2/' - % - 1 ) (2« - 1 ) j (59) 
^o(^(2/-%-l)(2«-l)) 



i' 

sin — / 2(-l 



(60) 
(61) 



In performing all these calculations, it is important to remem- 
ber that in current computing platforms, memory is both cheap 
and nicely managed so any number of the terms in these ar- 
rays can be calculated once or once per cycle and stored in 
advance. By employing distributed computing techniques, a 
dedicated CPU can spend it's time calculating these arrays 
while other parts of the program are running. Parallehzation 
of this algorithm will be discussed in detail below. 



E. Interaction Term 

The scattering term is written using the relaxation time ap- 
proximation 



Each processor is handling the lxH-2k problem, which im- 
plies a large number or matrix elements/equations. This num- 
ber is too large to satisfy the above conditions. In order to 
remedy this, we try to find a form of the matrix £2 that is block 
diagonal in k,-, similar to the ID problem, but with slightly 
different terms. We find that the potential term U cannot be 
expressed in such a way, and present an alternative method of 
solution, using approximation techniques, in section IVl 



A. Kinetic Matrix (Longitudinal) 

Following the procedure well defined in previous 
literature^ 01> we say [f{i,j)]i is written in vector 

r 1 ^ 

form as |^/li/l2/i.3 • ■■fi.NtJi,! ■ ■ ■ fij-ifijfij+i ■ ' J ^ 

so we can write that, for a given / the longitudinal drift term 
can be written in matrix form as (< and > denotes downwind 
and upwind differentiation, respectively.) 

fi.l.l 



2Tn 



^ + 1./ 



, (64) 



df{r,z,k,,kr) 



dt 



fo(r,z,k,,kr) 



t\S\K\dk>Jdk'Jo{ 



(62) 



where f(){r,z,K,k',.) is the equilibrium WDF. This term, dis- 
cretized, becomes 

\s-7] {n,ij,l) = P^"'''^''^ £ £ |(2/'-%-l)|/(«,,-,/,/')-i/(", 



and 



2T„ 



fi,\,l 



f. 



f. 



, (65) 



fo(nj.jj) 



I.%l!i',ti\i^l'-N,,-l)\Mn,iJ',n 



(63) 



IV. 2D MATRIX SETUP 

It is important that the matrix £2 = T + U + Sbe such that 
the limitations of modern computer systems are able to handle 
the above equations in a reasonable amount of time (days and 
weeks as opposed to months). We know of two general meth- 
ods that are able to efficiently solve the system of equations 
generated by the discretization of the WFE: Matrix inversion 
techniques and direct integration techniques (explicit Runge- 
Kutta-like, implicit BDF/ Adams, etc). For the matrix inver- 
sion methods, a necessary condition is that either the whole 
matrix can to be stored in RAM or can be broken up into parts 
that can be separately stored in RAM and solved for sequen- 
tially (i.e.: block diagonal matrix). The necessary condition 
for direct integration techniques is simply that the technique 
be stable and fast enough. 



being a diagonal square matrix of size N^, and C = ■ 

The values of = ±3, Tf' ~ =f4, = ±1 are defined by 
their position in the complete matrix, as stated by the second 
order differencing scheme and the boundaries. By denoting 
[/], as a vector of length A^^^ holding all the momentum (j) 
values of fij for a given i, the entire T ■ / term is written as 



T-/ =Cj 



( 



VL 



^1 ^0 

£7> C7> 

■^2 ^1 ^0 



tjT< 67< 
^0 ^1 



cjr< (jr< a7-< 
^0 ^\ ^2 

^0 ^1 

6r< 
^0 



rj7-> a7-> ej7-> 
^1 ^\ ^0 



[/] 



(66) 
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where Cj = {ij — Nk. — l)C, and T/ is a block tri-diagonal 
square matrix of rank N^Nt-^ . 

Concerning the radial terms, no matrix operations are 
needed. Each ring receives the Wigner function distribution 
from it's nearest inner and outer neighbor and the first order 
differencing scheme (described above) is used. 



/l, 1/1,2/1,3 • ■■f\.NtJ2,l ■ ■ -fij-lfi.jfi.j+l ■ ■ ■fN,.Nk- 

can write equation|^for a given j as 



, we 



B. Potential Matrix 

Now, when [/(/,y')]/' is expressed in vector form (for a given 
r and kr) as 



J 



y„,,(r,i-i,/') y,,,(/,i-2,/') 
y„,,(i,2-i,/') y,,,(/,2-2,/') 



V„j(Nk, - 1 - 1,/') V,j(i,Nt, - 1-2,/') 
V„.,(i,A'fc-l,/') y,,,(/,A'fc-2,/') 



Vnj{i,l-[Nk,-l\.l') V„j{i,l-Nk_,l') 

y„.,(i.2-[iVfc-i],/') y„,,{;,2-iv,^,/') 



[v(U')-aJ]„„ = c 

where V„,i{i,j-j',l') = f {n' ,l,l')^a(i' ,j ^ j')r{i,i\n,n') and C 

i' 

—V„j{i,j,l') and V„ i(i,0,l') = 0, the above becomes 

[v(u')-70] ,, = c 



V„j{i,Nt - 1 - [Nt, - 1]./') y,,,(i,iyt, - 1 -Nt^,l') 
V„j{iM - [Nk, - 1],/') V„j{iM, -NkJJ') 



/,.2 





V„.,(U,/') 



y„.;(;M,/') 




y„,,(i,%-2,/') y,,,(;,A't,-3,/') 

V^./C;,^*, -1,'') ^../(r,^*, -2,/') 
Note that V(/, /') is a A'^. x N^. anti-symmetric matrix. 



Since V„.i{iJ,l') °^ sin;, V„.i{i,~j,l') 



-V„,i{'M, -2,/') -V„j{i,Nt, - 1,?') 
-V„,;(/,Wfc,-3,/') -y„,,(;,Wt^-2,/') 







/-.2 



(67) 



By denoting [/], / as a vector of length A^^-.A^^.^ holding all the 
momentum (j) values of fiji for a given the entire U • / 



term is written as 



C. Interaction Matrix 



Whereas we previously wrote the discrete interaction term 
in equation|63] we rewrite it as 



S- / {n,i,j,l) = p{n,i)- -f{n,i,j,l) 

J T T 

Nk, 

p(«,/)= £ £ \{2l' -Nk,-l)\f{n,ij\l') 
/=l/'=l 

fo{n,iJ,l) 



P{n,iJ,k) = 



i/;ii!^tii(2''-A^-t.-i)i/o(«,',/,n 



(69) 



V(l,l) V(l,2) 
V(2,l) V(2,2) 



V(l,%) 
V(2,A^,J 



so that U is a square matrix of rank N^Nk.Nkr- 





r [/]i,i 1 




[/]l.2 






(68) 



where p (n, /) is the density of the previous time step. 

By writing/(/,y) as a vector, 

[/i, 1/1,2/1,3 • --fLNkfiA ■ ■■fi.j-ifijfij+i ■ ■ ■/a'v.a'J ^ and 
by denoting [/], as a vector of length holding all the 

momentum (j) values of fij for a given i, the entire S ■ / term 
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IS written as 



n,l 



s-f 



S{1) 
5(2) 





_ 1 

t,l X 

■ 

■ 



[PW]2 

[P]ivjj3]iv, 









[fh 




- [f]N, - 



(70) 



where S is a square block diagonal matrix of rank NzNk- and 
S is a N^N^^ vector. 

D. Boundary Conditions 

By denoting [/],■ as a vector of length A'^^ holding all the 
momentum (j) values of the Fermi distribution fFermiU) for a 
given i, the longitudinal boundary equations (l31ltol?^ can be 
written as 



B = Cj 



B^fh 



Bt [An,-, 
Bf [An, 



(71) 



Bn" [f]i being a vector of size A^^.,, B a vector of size N^N^., and 

Cj = ^li {y-Nk, - 1). The values ofsf = ±2, flf = t1 
are defined by equations |^ to |53 as stated by the second 
order differencing scheme at the boundaries. 



V. METHODS OF SOLUTION 

A. Parallelization 

So far, we have dealt almost exclusively with the ID space, 
2D momentum transport problem (lx+2k). As stated above, 
we are working under the assumption that longitudinal trans- 
port is more dominant than radial transport. This allows the 
total transport to be calculated in two steps: (1) transport the 
particles in the longitudinal direction in each shell separately 
(lx+2k), then (2) each shell exchanges particles with its near- 
est neighbor. This latter step is where we employ parallel pro- 
cessing techniques. The lxH-2k problem is performed on P 
processors, where P is the number of cylindrical shells into 
which we have divided up the RTD. Once each shell has ad- 
vanced a given amount of time, then communication between 
shells (radial drift) can commence. As described in Section 
nil CI a central differencing scheme (CDS) is used and shown 
in eauations l4 11431 The boundaries consist of the material ex- 
ternal to the shell and the innermost shell. As per equations 
144 1471 and explained in section IIII CI the exterior boundary 
defines the device. For example, if the device is a mesa RTD, 



and there is nothing but vacuum outside the shell, one should 
choose a boundary shell that injects into the outermost shell 
the same particles that the outermost shell ejected (keeping 
the momenta the same). If the device is a slab with a circular 
contact for an emitter, then the material outside the cylinder is 
in equilibrium. The boundary shell would be chosen to reflect 
this. 



B. Potential Transform 

As explained is section IIVI and seen in equation |68j the 
Wigner integral (potential term) is the only term that can not 
be made diagonal in kr, leading to a full matrix that is too 
big to store in memory. If one uses direct integration methods 
(01 )' the number of terms becomes large enough to make the 
problem intractable. We now describe a method of using the 
Fourier transform property of the Wigner integral to eliminate 
the kr dependence inherent in eauationFTTI 

We need to have either the matrix i2 = T^^ + U + S be able to 
be stored in the RAM of present day computers (for the matrix 
methods), or limit the number of simultaneous equations to 
solve ( the direct integration methods). The idea behind our 
method of solution of the 2D transport equation. 



dt 



ilf B, 



(72) 



is that if the matrix Q. is block diagonal in kr then one can 
progress through all the values of kr solving the matrix equa- 
tion at fixed values of kr each time, effectively reducing the 
simulation to a series of ID problems. Unfortunately, the 
equations for some of the matrix operators are not block 
diagonal in kr in their present form. Some manipulation 
will be needed to obtain block diagonal terms. The drift 
term is already independent of kr, and the scattering term is 
coupled to kr via off-diagonal elements due to the integral 
f\k'r\'^K I dk'^f{r,z,k'.^k'r). This integral is just the IDdensity 
in the shell Pr{z), which can be calculated at the beginning of 
the time step. This approximation makes the scattering term 
into the desired diagonal matrix without losing too much de- 
tail. The potential term, however, is not so simple. In equation 
( BSt . the term ^(r,z,p, Q, as defined in equation (tTTTi makes 
the matrix operator U a full matrix. We will now outline a 
method to circumvent this problem below. 

In order to solve the transport equation, we use an implicit 
method (first proposed in QIl) by rewriting equation 1721 as 
(dropping the z index from the potential matrix) 



f f f+f 

^ = (T + U + S)^-B, 



(73) 



where f means the new (next) value of f in time. Rewriting it 

as 



1-^(T + U + S) 



(f+f) =2f+BAf, (74) 



and making the approximation (accepting error in terms of 
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order A/^) 

l-y(T + U + S): 

allows us to write that 



l-y(T + S) 



1-|(T + S) 



• (f+f) =2f + BAf. 



For convenience, we will rewrite this as 

at = Ooaut = b 

where we have defined 
b 



T = 



f 



2(f + BT) 
Af 
2" 

(1-tU) 
(1-T[T + S]) 
(f+f), 



(75) 

(76) 

(77) 

(78) 
(79) 

(80) 
(81) 
(82) 



The solution to this equation involves solving two matrix 
equations. 

1 . Solve iioF = b for F (quick) 

2. Solve £2{/f' = F for f (impractical) 

The matrix i2[/ is still too big to store in memory and, 
consequently, solving i2j/f = F for f is not practical for 
modern computers. By making the approximation that led to 
equation j76> . we are able to solve i2f/f = F separately, 
allowing us to reformulate it in a form more suitable for 
computation. 

Recall that the Wigner integral, U • f , was derived by taking 
the Fourier transforms of the Greens function, \ 8, sec 6.4]. 
From the last term in equation^! can be written out in full as 



dz e 



f„-ilk-j. 



dk 



t +i2kU' 



''J \r'\dr'Jo{2krr')f'{z,z',r,r')x 
k'AdklJo(2k'/)f{r,z,k[,k'^) 



FrF^'f, 



(83) 



rlTt 

'r{z,^,r,r') = rfe {y(z + z',r + r' cos e)-y(z-z',r-z' cose)}. 



(84) 

The Fourier Bessel transform,F, and its inverse, F^', are 
given by 

¥ = ¥{kr,k,\r ,z) = I dze-'-''^- j dr'\r'\jo (2^^/) (85) 

F-^ = F-\kr,k,-r'.z') = — '-^ / / dk,\kr\Jo{2ky). 
{2nY J J 

(86) 

When we define 

^(x,y) = / d^e^'^^y f{^M = F-'f. (87) 

we can rewrite i2{/f = F as 



(1-TFrF-i) (Fg') = F, 



(88) 



g being the Fourier Transform of f and g' = g + g. Some ma- 
nipulation gives 



where 
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F^'F. 



(89) 



(90) 
(91) 



The new procedure is to solve for the new equation, i2[,g' = 7. 

The term ^(r,z,p, ^), as defined in equation ( fill , can now 
be written as 

= y t/3k'e2ikV(r,z,k) 

= (2;r)3g(r,z,p,C,0). (92) 
This lets us define H'^ = [l — yU'] in terms of 
1 



U'g: 



dpdQp\sm{2k,QJo{2krp) x 



y{r,z,p,Q{2n)'g{r,z,p,C,0), 



(93) 



or, in discrete form, as (completing equation ( I57t ) 

[U'g] («,/,;,^) =+^Ar2Az£^(«,/,n',/',^,j)) 



n' ,r 



g{n,i,n ,i ,m ) 
The WDF is recovered by 



(94) 



f{r,z,k,-,L) = \kr\dX(i,f{r,z,kr,k-,X(l>) 



2n 



kr\dX^ d'ye^'''yg{r,z,y) 



k,-\ J dp\p\Ja{2krP) I dCe 



-'2ik-i; 



(95) 



The complete process to solve the WFE equation is given in 
the following steps 

1. Solve HoF = b for F 

2. Define 7= F^'F. 

3. Solve iiyg' = 7 for g' 

4. Recover f from f = Fg' 
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C. Aim and shoot 

We find that using direct integration for a system of linear 
equations grew slower and less stable as the number of simul- 
taneous equations grew. For example, a Ifs time step involved 
hundreds of iterations were involved, each one taking about 90 
seconds. When using an implicit matrix method, as is standard 
in the ID codes, the matrices were too big to store. Because 
we have now, using the approximation in equation ^5] sepa- 
rated the operator Q. into the product of a drift/scattering oper- 
ator, Q.0, and a potential operator, £2[/, a combination of these 
two methods can be used. As we shown, by itself, CIq can be 
block diagonal and sparse enough to be easily solved by ei- 
ther of the two methods. Also, by itself, £2j/ can be rewritten 
in a form that also can be easily solved by either of the two 
methods. 

The following method allows us to take advantage of this 
split is a way that is analogous to the accelerated convergence 
method used to obtain steady state solutions in the ID problem 
We have dubbed this a "aim and shoot" approach. The 
static potential term is calculated on a time scale of 5t = 0.01 
fs for only one step, then held constant while the system 
evolves for At = Ifs. This is the aim part. The shoot part 
involves the drift and scattering matrices being solved implic- 
itly by matrix inversion, with the potential term static. This is 
fine for a steady-state solution, but it is still uncertain if this 
method will give the proper transient behavior of the system. 




Figure 2: Initial WDF as Boundary Value, (a) Longitudinal phase 
space for one fc,. slice, (b) All k,- slices put together. 




position 



longitudinal momentum 



Figure 3: Initial WDF as a fanciful test distribution. 



VI. IMPLEMENTATION AND SIMULATION RESULTS 

In the preceding sections a computational method of solv- 
ing for both the time dependent and steady state two dimen- 
sional Wigner function transport equation was presented. The 
2D equations and computational method were derived for the 
case of longitudinal transport through a cylinder while taking 
account of the effects of radial momentum in addition to the 
longitudinal momentum. As previously stated, the numerical 
solution is broken into two parts: (1) transport in the longitu- 
dinal direction then (2) transport in the radial direction. The 
cylinder is divided in to a number of concentric cylindrical 
shells in which the longitudinal transport takes place as in the 
ID problem, but with the inclusion of radial momentum. The 
radial transport involves a simple exchange of particles (dic- 
tated by the newly calculated radial momentum). Sine this 
step is computationally trivial, this work was concerned with 
the former step: A ID space and 2D momentum transport 
problem (lxH-2k). Below we present some proof-of-principle 
simulation results obtained using the methods developed in 
section |V] From these results, the future usefulness of each 
of these methods, in light of cuiTent computing trends, will be 
discussed. 

This simulations were caiTied mostly on Linux worksta- 
tions (2GHz Pentium4, 1.7GHz Pentium4s and 1.2GHz AMD 
AthlonMPs). For phase space, our solution method is most 
easily formulated when = Nk. and Nr = Nk, since the 
Fourier transforms between momentum space and displace- 
ment space must be on the same lattice. Simulations were per- 



formed on longitudinal phase space grid sizes of A^, = A'^^. = 
96 and A^, = A'^. = 48 with radial phase space grid sizes of 
Nr = Nk, = 2,4,8, 16,20. In all of the simulations presented, 
the device is constructed of bulk n-doped GaAs with a RTS of 
undoped GaAs /Alo^Gao jAs with a baiTier potential is 0.3 eV. 
The device temperature is 77K, the electron effective mass is 
0.0667mo and The donor density is 2 x 10^*^ cm^^. The cylin- 
der has dimensions of 1000 A in both length (z) and diameter 
(r). For most of the simulations the active RTS region is ap- 
proximately 180A. The well, spacer and barrier lengths will 
be given for each simulation. 

Each numerical experiment has been carried out in the fol- 
lowing way. First, the cylinder is populated with electrons 
according to one of three specific WDF: (a) f{z,kz,kr) = 
0, corresponding to no excess electrons in the cylinder, (b) 
f{z,kz,kr) — f Fermi, Corresponding to the same distribution 
of the metallic leads (figure|2jl, and (c) a fanciful distribution 
coiTesponding to electrons mostly at the center of phase space 
(figure 13. To be exact, an error was made in preparing for 
case (b). We meant to set the WDF each kr slice to the Fermi 
distribution of the boundary of that specific kr slice. This way 
the integral of the total WDF would be unity, as expected. In- 
stead, we accidentally normalized the WDF of each A:rSlice 
such that the integral of the WDF in each slice in unity. We 
have kept this eiTor here for the reason that the simulations 
illustrates that the system will adjust itself and still tend to- 
ward the expected result. Next, at zero bias, the system 
is allowed to evolve with scattering turned off. The system 
is allowed to evolve for a suitable time, until it settles into a 
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steady state, and then scattering is turned on. During this time, 
the previous step's value of the WDF is used as the "equilib- 
rium value" needed in the relaxation time approximation of 
the present step. Once again, the system is allowed to evolve 
for a suitable time, and at this time, the current WDF is set to 
be the "equilibrium value" for the rest of the evolution, which 
continues until the system settles into a steady state. By ob- 
serving this time evolution of each of the three cases, we can 
determine how well a given method behaves as well as gather 
timing and other information. 

One important item of note is that the examples below are 
performed with A^^^ = Ni^, ~ 96 and A^^ = N^^ = 2. In effect we 
are including only one positive and one negative radial mo- 
mentum value. While this is fine for testing purposes where 
we are basically reducing the simulation to a ID problem, for 
a real 2D problem A'^^ and Nr should be large enough (-16) 
to encompass a phase space greater than the radial Fermi mo- 
mentum of the material outside the cylinder The computa- 
tional issue is that while A^^ = — 2 can be solved in under 
20min for a 2000fs run at time steps of Ifs, the addition of 
more radial points increases the time dramatically (this will 
be mentioned below). As a result, the phase space plots given 
are showing only one value of the radial momentum (the pos- 
itive momentum) since they are symmetric in this case. 

In order to follow the time evolution of this system, we 
have employed two different types of integration methods: 
Direct integration (explicit) and Matrix solvers (implicit). Di- 
rect integration is done using a pre-packaged integrator called 
ROCK4|7], which is a 4''' order Runge-Kutta like integrator 
for a system of equations. With ROCK4, there is no need to 
compute and store the right hand side of the discretized equa- 
tion as a general band matrix, and consequently, no need for 
a matrix inversion of the time evolution operator which have 
been used in previous simulations. Instead, all that is needed 
is to calculate the right hand side on the fly. Implicit integra- 
tion is done by using LAPACK[9J to solve the matrices. 





Potential 


(/ 7^0 


U 


= 




# time steps 


200Ar 


2000A/ 


200Af 


2000A/ 


Initial 


Zero 


0m7.428s 


Om5 1.376s 


lm38.178s 


16m37.437s 


WDF 


PseudoFermi 


0m7.349s 


Om5 1.762s 


lm36.676s 


15m45.484s 


Value 


Central 


0m7.396s 


Om52.361s 


lm35.993s 


15m38.711s 



Table I: Timings for the Matrix Split runs 

of equation ^5] and split the matrix Q. into a product of a 
drift/scattering matrix and a potential matrix (equation I76>. 
This can be written as 



= £2oi2t/f' = b, 



(96) 



The solution to this equation involves solving two matrix 
equations, one for the drift/scattering 



iinr = b, 



and one for the potential 



(97) 



(98) 



In solving the potential equation, we must perform an inverse 
ti-ansform, 7 = F^'P, and then a transform, f = Fg'. 

Figures l4l through 1151 illustrates the time evolution of our 
device up to 2000fs. Each figure is a snapshot in the evolu- 
tion of three WDFs whose initial values are one of the three 
discussed above. They will be referred to the Zero, Pseud- 
oFermi and Central initial values. In each figure, the left col- 
umn shows the WDF plotted in 3D, while the right column 
shows a contour plot of the WDF. We will examine two dis- 
tinct sets of cases to illustrate our method. The first will have 
the potential term set to zero, corresponding to no coulombic 
interaction. The second will include the coulombic effects. 



A. Direct Integration Methods 

Not much will be said for the explicit method since, on the 
same computer as the runs below, after 36 minutes the method 
progressed only to a time of 14fs. When the number of ra- 
dial grid points is increased from 2 up to 8, the number of 
equations increases 256 times the original number. This fact 
renders ROCK4 useless for any future 2D simulations. Re- 
cently, we have been introduced] 10| to implicit direct integra- 
tion methods (BDF/ Adams) and Newton solvers that, so far, 
outperform the ROCK4 method for the ID simulations. We 
have yet to include this method in our 2D simulations. The 
next phase of the ongoing research is to see how such meth- 
ods compare to what we will present below. 



B. Matrix Split 

The so-called matrix split refers to the method of section 
IV Bl which is when take the WFE, apply the approximation 



1. Potential "Turned Off" 

Figures|3through|9lshow the evolution of the electron dis- 
tribution in our device without any coulomb interactions. We 
see that for each of the three WDFs, the system behaves as 
expected. With a zero initial WDF, we see the electrons move 
in from either end of the device with the high momentum car- 
riers further in than the low momentum carriers. As time pro- 
gresses, the phase space fills with carriers and the WDF tends 
towards the distribution of the boundaries, namely, the Fermi 
distribution. The next case, that of the initial WDF set to the 
PseudoFermi distribution, shows the excess carriers leaving 
the system. Ultimately, the phase space once again tends to- 
wards the distribution of the boundaries. In the final case, the 
Central distribution, we see the carriers moving in from the 
boundary, as in the first case. At the same time, the central 
distribution itself relaxes, with the positive momentum carri- 
ers moving one direction and the negative momentum carri- 
ers the other way. Eventually, the central distribution relaxes 
while the boundary carriers move in. Once again, the end re- 
sult tends towards the boundary distribution. 
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All three of the different initial WDFs evolve towards 
the same, expected, final WDF. This simply shows that the 
drift/scattering part of the simulation works as expected, 
which is expected. What we can learn from this exercise is 
the CPU time to calculate the drift and scattering terms up to 
200 and 2000 fs (Table|D. As stated above, this is the first of 
two matrix equations that must be solved. We see from the ta- 
ble that this is not where most of the CPU will spend its time. 
Rather, the potential term will take the bulk of the computing 
time. Next, we see how the simulation behaves when this term 
is turned on. 

DUE TO ARXIV SIZE CONSIDERATIONS THESE GRAPHS 
ARE PLACED IN A SEPARATE FILE AVAILABLE FOR 
DOWNLOAD ALONG WITH THIS 



Figure 4: Neglecting coulomb interactions: Surface and Contour 
plots of the WDF at=10fs for the Initial WDF of (a) zero, (b) Pseud- 
oFermi, (c) Central 
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Figure 5: Neglecting coulomb interactions: Surface and Contour 
plots of the WDF at=50fs for the Initial WDF of (a) zero, (b) Pseud- 
oFermi, (c) Central 
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Figure 6: Neglecting coulomb interactions: Surface and Contour 
plots of the WDF at=100fs for the Initial WDF of (a) zero, (b) Pseud- 
oFermi, (c) Central 
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Figure 7: Neglecting coulomb interactions: Surface and Contour 
plots of the WDF at=200fs for the Initial WDF of (a) zero, (b) Pseud- 
oFermi, (c) Central 
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Figure 8: Neglecting coulomb interactions: Surface and Contour 
plots of the WDF at=1000fs for the Initial WDF of (a) zero, (b) 
PseudoFermi, (c) Central 
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Figure 9: Neglecting coulomb interactions: Surface and Contour 
plots of the WDF at=2000fs for the Initial WDF of (a) zero, (b) 
PseudoFermi, (c) Central 



2. Potential "Turned On" 



Figures [TUI through 1151 show the evolution of the electron 
distribution in our device including coulomb interactions. We 
can compare the simulations of these three cases with coulom- 
bic interactions to those above, where coulombic interactions 
were ignored. By following the evolution of the first case 
(zero initial WDF), we see how the carriers interact with the 
barriers as they move towards the center of the device. We 
also see (more clearly in the contour plots) how the carriers 
interact with each other, spreading out slightly as the progress 
inwards. In figures [21 and [O] we begin to see the interaction 
of the reflected carriers with the incoming carriers. This effect 
grows as the system evolves, which is evident in the dark/light 
patterns along the momentum axis. The same effects are seen 
in the other two cases, with the exception that they begin with 
the unlikely distribution having carriers in the barriers. We 
see these carriers begin ejected from the barrier region at high 
momentum at first, then at later times these cases evolve to 
same result as the first case, where we see the expected WDF. 

As we increase A^^ — N^^ from 2 to 4, a 200 fs run increases 
from about lm40s to 5min (a factor of 3). Projecting this to 
2000fs, we see the simulation would take approximately 47 
min to complete. So far, this is not unreasonable, since for 
a given bias point, a 2000fs run is enough to assure conver- 
gence. A 90 point IV curve (including the reverse sweep), 
would take about 70 hours (about 3 days). Increase Nr = Nt, 
to 8 and a 200fs run takes 20min ( 4 x A^^ = 4, or 12 x A^^ = 2). 
A 2000fs run will take 3|hours, which means a 90 point IV 
curve will take 12.5 days. A beginning run at Nr ~ N^^ — 16 
returns a time rate of 38 sec/fs. At that rate, a 2000fs run 
would take 21.1 hours and a 90 point IV curve almost 80 days. 
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Figure 10: Including coulomb interactions: Surface and Contour 
plots of the WDF at=10fs for the Initial WDF of (a) zero, (b) Pseud- 
oFermi, (c) Central 
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Figure 11: Including coulomb interactions: Surface and Contour 
plots of the WDF at=50fs for the Initial WDF of (a) zero, (b) Pseud- 
oFermi, (c) Central 
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Figure 12: Including coulomb interactions: Surface and Contour 
plots of the WDF at=100fs for the Initial WDF of (a) zero, (b) Pseud- 
oFermi, (c) Central 
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Figure 13: Including coulomb interactions: Surface and Contour 
plots of the WDF at=200fs for the Initial WDF of (a) zero, (b) Pseud- 
oFermi, (c) Central 
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Figure 14: Including coulomb interactions: Surface and Contour 
plots of the WDF at=1000fs for the Initial WDF of (a) zero, (b) 
PseudoFermi, (c) Central 
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lions have been restricted to ID drift/scattering and ID poten- 
tial. Any comparison to a real device must answer the ques- 
tion of how important both radial drift/scattering and a non 
ID Coulombic charge density are. We are now prepared to 
examine these questions. 

The computational hurdles of solving a 2D WFE have been 
identified: (1) Not all the matrices are sparse enough to fit into 
the RAM of present day computers, and (2) Direct integration 
becomes more time consuming as the number of simultaneous 
equations to solve grows large. Some solutions of these hur- 
dles have been described and a workable way of numerically 
simulating a RTS that exhibits cylindrical symmetry has been 
given. If one can separate the operator i2 = Tj + U + S into 
the product of a drift/scattering operator and a potential op- 
erator, then each matrix can be made sparse (overcoming the 
first hurdle) and/or block diagonal (the second hurdle). This 
also allows the separation of the transport problem into a two 
step problem, drift/scatter and potential computations. An ad- 
ditional time saving method, based on this split, is introduced. 
Finally, but taking advantage of the Fourier transform nature 
of the potential term, one can decrease the operator size sub- 
stantially further. 

Time evolution simulations based on these method were 
then presented for three different cases. Each case lead to nu- 
merical results consistent with expectations. To the author's 
knowledge, this is the first proof-of-principle of a ID space, 
2D momentum simulation. At the present time, a full tran- 
sient treatment of a forward/backwards bias sweep would take 
upwards of 3 months. The work shown here still must be nu- 
merically scrutinized in order to shorten the computer times 
involved. We are beginning to work with a group of numerical 
experts to enhance the PDE solvers shown here. One impor- 
tant item to remember is that the computational hardware is 
still following Moore's law, that is, approximately doubling in 
speed every 18 months. In a few years, the techniques shown 
here, which push current technology to the limit, will prove to 
be even more feasible in the foreseeable future. 

It is our belief that the methods outhned in this paper will 
finally allow a fuU treatment of the RTS transport problem. 



Figure 15: Including coulomb interactions: Surface and Contour 
plots of the WDF at=2000fs for the Initial WDF of (a) zero, (b) 
PseudoFermi, (c) Central 



VII. CONCLUSIONS 

We have shown here a a 2D drift/scattering and 3D potential 
form of the Wigner function transport equation for the case of 
a cylindrical device. This is an important step in Wigner func- 
tion simulations of electronic transport since previous simula- 
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